library(psych)
library(nlme)
library(effects)
library(ggplot2)
library(metafor)
library(grid)

############################################################################################################################################
############################################################# MTurk 2016 Study #############################################################
############################################################################################################################################

setwd("C:/Users/SCote/Dropbox/Research/Projects/Social Class by Inequality Effects/PNAS letter")
data <- read.csv("MTurk.2016.Study.csv")

with(data, describe(cbind(personalincome,householdincome,stateGini,DGgenerosity)))

grpmeans<-aggregate(cbind(data$householdincome,data$personalincome),list(data$stateliveFIPS),mean,na.rm=TRUE)
names(grpmeans)<-c("stateliveFIPS","state.mean.householdincome","state.mean.personalincome")
dataHLM<-merge(data,grpmeans,by="stateliveFIPS")
dataHLM$householdincome.groupc<-dataHLM$householdincome-dataHLM$state.mean.householdincome
dataHLM$personalincome.groupc<-dataHLM$personalincome-dataHLM$state.mean.personalincome
dataHLM$householdincome.groupc.div <- dataHLM$householdincome.groupc/10000
dataHLM$personalincome.groupc.div <- dataHLM$personalincome.groupc/10000

# analyses with household income
summary(HLMmodel.household <- lme(DGgenerosity ~ householdincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
with(dataHLM,describe(cbind(householdincome.groupc.div,c.stateGini,DGgenerosity)))
mean(dataHLM$c.stateGini); sd(dataHLM$c.stateGini); max(dataHLM$c.stateGini); min(dataHLM$c.stateGini)
donationSESGiniHLM1<- data.frame(effect("householdincome.groupc.div*c.stateGini", HLMmodel.household, xlevels=list(c.stateGini=c(-0.04006078,0.009688921,0.04653922))))
ggplot(donationSESGiniHLM1) + geom_line(aes(householdincome.groupc.div, fit, color = c.stateGini, linetype=factor(c.stateGini)), size= 1.25, )+ scale_colour_gradient(low = "blue", high = "red") + theme_bw()+
  coord_cartesian(xlim=c(-10,50), ylim=c(0,10)) +
  scale_linetype_manual(values=c("dotted", "dashed", "solid"))+
  theme(legend.position='none', panel.grid.major = element_blank(),
        text = element_text(size=16),
        panel.grid.minor = element_blank(),)+ xlab("Household Income")+ylab("Dictator Game Giving")
# verify that results are comparable with grand-mean centered income
dataHLM$householdincome.div<-dataHLM$householdincome/10000
summary(HLMmodel.household.div <- lme(DGgenerosity ~ householdincome.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

# analyses with personal income
summary(HLMmodel.personal <- lme(DGgenerosity ~ personalincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
with(dataHLM,describe(cbind(personalincome.groupc.div,c.stateGini,DGgenerosity)))
mean(dataHLM$c.stateGini); sd(dataHLM$c.stateGini); max(dataHLM$c.stateGini); min(dataHLM$c.stateGini)
donationSESGiniHLM2<- data.frame(effect("personalincome.groupc.div*c.stateGini", HLMmodel.personal, xlevels=list(c.stateGini=c(-0.04006078,0.009688921,0.04653922))))
ggplot(donationSESGiniHLM2) + geom_line(aes(personalincome.groupc.div, fit, color = c.stateGini, linetype=factor(c.stateGini)), size= 1.25, )+ scale_colour_gradient(low = "blue", high = "red") + theme_bw()+
  coord_cartesian(xlim=c(-10,15), ylim=c(0,10)) +
  scale_linetype_manual(values=c("dotted", "dashed", "solid"))+
  theme(legend.position='none', panel.grid.major = element_blank(),
        text = element_text(size=16),
        panel.grid.minor = element_blank(),)+ xlab("Personal Income")+ylab("Dictator Game Giving")
# verify that results are comparable with grand-mean centered income
dataHLM$personalincome.div<-dataHLM$personalincome/10000
summary(HLMmodel.personal.div <- lme(DGgenerosity ~ personalincome.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

############################################################################################################################################
####################################################### Human Cooperation Lab Study  #######################################################
############################################################################################################################################

data <- read.csv("Human.Cooperation.Lab.Study.csv")

grpmeans<-aggregate(data$personalincome,list(data$stateliveFIPS),mean,na.rm=TRUE)
names(grpmeans)<-c("stateliveFIPS","state.mean.personalincome")
dataHLM<-merge(data,grpmeans,by="stateliveFIPS")
dataHLM$personalincome.groupc<-dataHLM$personalincome-dataHLM$state.mean.personalincome
dataHLM$personalincome.groupc.div <- dataHLM$personalincome.groupc/10000

# analyses with personal income
summary(HLMmodel.personal <- lme(DGgenerosity ~ personalincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
with(dataHLM,describe(cbind(personalincome.groupc.div,c.stateGini,DGgenerosity)))
mean(dataHLM$c.stateGini); sd(dataHLM$c.stateGini); max(dataHLM$c.stateGini); min(dataHLM$c.stateGini)
donationSESGiniHLM2<- data.frame(effect("personalincome.groupc.div*c.stateGini", HLMmodel.personal, xlevels=list(c.stateGini=c(-0.05826078,0.00909524,0.07563922))))
ggplot(donationSESGiniHLM2) + geom_line(aes(personalincome.groupc.div, fit, color = c.stateGini, linetype=factor(c.stateGini)), size= 1.25, )+ scale_colour_gradient(low = "blue", high = "red") + theme_bw()+
  coord_cartesian(xlim=c(-10,25), ylim=c(0,10)) +
  scale_linetype_manual(values=c("dotted", "dashed", "solid"))+
  theme(legend.position='none', panel.grid.major = element_blank(),
        text = element_text(size=16),
        panel.grid.minor = element_blank(),)+ xlab("Personal Income")+ylab("Dictator Game Giving")
# verify that results are comparable with grand-mean centered income
dataHLM$personalincome.div<-dataHLM$personalincome/10000
summary(HLMmodel.personal.div <- lme(DGgenerosity ~ personalincome.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

############################################################################################################################################
##################################################### Pre-Registered MTurk 2019 Study  #####################################################
############################################################################################################################################

data <- read.csv("Preregistered.MTurk.2019.Study.csv")

# correlations between measures of generosity
with(data, corr.test(cbind(DGgenerosity,volunteering,charitydonations)))
with(data, cor.test(DGgenerosity,volunteering))

grpmeans<-aggregate(cbind(data$householdincome,data$personalincome),list(data$stateliveFIPS),mean,na.rm=TRUE)
names(grpmeans)<-c("stateliveFIPS","state.mean.householdincome","state.mean.personalincome")
dataHLM<-merge(data,grpmeans,by="stateliveFIPS")
dataHLM$householdincome.groupc<-dataHLM$householdincome-dataHLM$state.mean.householdincome
dataHLM$personalincome.groupc<-dataHLM$personalincome-dataHLM$state.mean.personalincome
dataHLM$householdincome.groupc.div <- dataHLM$householdincome.groupc/10000
dataHLM$personalincome.groupc.div <- dataHLM$personalincome.groupc/10000

# analyses with household income
summary(HLMmodel.household <- lme(DGgenerosity ~ householdincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
with(dataHLM,describe(cbind(householdincome.groupc.div,c.stateGini,DGgenerosity)))
mean(dataHLM$c.stateGini); sd(dataHLM$c.stateGini); max(dataHLM$c.stateGini); min(dataHLM$c.stateGini)
donationSESGiniHLM1<- data.frame(effect("householdincome.groupc.div*c.stateGini", HLMmodel.household, xlevels=list(c.stateGini=c(-0.05826078,0.008666553,0.07563922))))
ggplot(donationSESGiniHLM1) + geom_line(aes(householdincome.groupc.div, fit, color = c.stateGini, linetype=factor(c.stateGini)), size= 1.25, )+ scale_colour_gradient(low = "blue", high = "red") + theme_bw()+
  coord_cartesian(xlim=c(-10,50), ylim=c(0,10)) +
  scale_linetype_manual(values=c("dotted", "dashed", "solid"))+
  theme(legend.position='none', panel.grid.major = element_blank(),
        text = element_text(size=16),
        panel.grid.minor = element_blank(),)+ xlab("Household Income")+ylab("Dictator Game Giving")
# verify that results are comparable with grand-mean centered income
dataHLM$householdincome.div<-dataHLM$householdincome/10000
summary(HLMmodel.household.div <- lme(DGgenerosity ~ householdincome.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

# analyses with personal income
summary(HLMmodel.personal <- lme(DGgenerosity ~ personalincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
with(dataHLM,describe(cbind(personalincome.groupc.div,c.stateGini,DGgenerosity)))
mean(dataHLM$c.stateGini); sd(dataHLM$c.stateGini); max(dataHLM$c.stateGini); min(dataHLM$c.stateGini)
donationSESGiniHLM2<- data.frame(effect("personalincome.groupc.div*c.stateGini", HLMmodel.personal, xlevels=list(c.stateGini=c(-0.05826078,0.008662646,0.07563922))))
ggplot(donationSESGiniHLM2) + geom_line(aes(personalincome.groupc.div, fit, color = c.stateGini, linetype=factor(c.stateGini)), size= 1.25, )+ scale_colour_gradient(low = "blue", high = "red") + theme_bw()+
  coord_cartesian(xlim=c(-10,50), ylim=c(0,10)) +
  scale_linetype_manual(values=c("dotted", "dashed", "solid"))+
  theme(legend.position='none', panel.grid.major = element_blank(),
        text = element_text(size=16),
        panel.grid.minor = element_blank(),)+ xlab("Personal Income")+ylab("Dictator Game Giving")
# verify that results are comparable with grand-mean centered income
dataHLM$personalincome.div<-dataHLM$personalincome/10000
summary(HLMmodel.personal.div <- lme(DGgenerosity ~ personalincome.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

# analyses with charitable donations
summary(HLMmodel.household.charitydonations <- lme(charitydonations ~ householdincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
summary(HLMmodel.personal.charitydonations <- lme(charitydonations ~ personalincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

# analyses with volunteering
summary(HLMmodel.household.volunteering <- lme(volunteering ~ householdincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))
summary(HLMmodel.personal.volunteering <- lme(volunteering ~ personalincome.groupc.div*c.stateGini, random=~1|stateliveFIPS, na.action="na.exclude", data=dataHLM))

############################################################################################################################################
############################################################## Meta Analysis  ##############################################################
############################################################################################################################################

#### meta-analysis of household income

Study <- c(1,2,3)
Label <- c('Original PNAS 2015 Study','MTurk 2016 Study','Pre-Registered MTurk 2019 Study')
SampleSize<- c('N=1498','N=507','N=3806')
meta.data <-data.frame(Study,Label,SampleSize)

Interaction<-c(-1.56,-2.49,-.23) # coefficients
Error<-c(.54,1.52,.27) # standard errors
fixed.meta1 = rma(yi = Interaction, sei = Error, method = "FE")
summary(fixed.meta1)
forest(fixed.meta1, slab = paste(meta.data$Label, meta.data$SampleSize, sep = ", "))
grid.text("Meta-Analysis of Interaction Coefficients", .5, .9, gp=gpar(fontsize=8,cex=2))
grid.text("Across Three Studies", .5, .85, gp=gpar(fontsize=8,cex=2))

#### meta-analysis of personal income

Study <- c(1,2,3)
Label <- c('MTurk 2016 Study','Human Cooperation Laboratory Study','Pre-Registered MTurk 2019 Study')
SampleSize<- c('N=507','N=3696','N=3802')
meta.data <-data.frame(Study,Label,SampleSize)

Interaction<-c(-5.55,-.91,.22) # coefficients
Error<-c(2.34,.45,.41) # standard errors
fixed.meta1 = rma(yi = Interaction, sei = Error, method = "FE")
summary(fixed.meta1)
forest(fixed.meta1, slab = paste(meta.data$Label, meta.data$SampleSize, sep = ", "))
grid.text("Meta-Analysis of Interaction Coefficients", .5, .9, gp=gpar(fontsize=8,cex=2))
grid.text("Across Three Studies", .5, .85, gp=gpar(fontsize=8,cex=2))

